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ABSTRACT 


Military electronics are frequently operated in excessively confined spaces aboard 





ships and aircraft. This limited space impacts significantly on the space available 
for cooling equipment. The optimal solution is the development of one universal, 
modular rack for shipboard or aviation use. With a modular design, upgrades to 
equipment could also be accompanied by an upgrade to the cooling rack itself with 
very little additional cost or difficulty. A water cooled avionics rack can provide suf- 
| ficient cooling for any piece or combination of avionics equipment if enough water 
flow paths are used, the water is at the appropriated temperature and the water 
is properly distributed within the passages. To determine if the cooling medium, 
water, is sufficiently distributed within a modular cooling rack, an analysis of the 
flow and pressure distribution of the coolant is required. This thesis presents a 
computer code that has been developed as an initial step in the total design of a 
modular cooling rack for avionics equipment. In itself, the code details a specific 
design technique and allows for the determination of whether a proposed configu- 
ration, including source location, characteristics of the cooling water, and the size 
and shape of the proposed flow passages will indeed provided a proper distribution 


of the coolant. 
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I. INTRODUCTION 
A. GENERAL 


This thesis describes the development of a computer code to analyze the flow 
and pressure distribution of water in an avionics cooling rack. This code allows for 
variable rack design and flow sources. The code is written in FORTRAN 77. In 


addition, the code can be run on any IBM or IBM compatible personal computer. 


B. BACKGROUND 


Extremely high temperatures are a primary cause of avionics equipment un- 
reliability. The origin of the thermal problem is in the continuous effort to develop 
lighter and smaller components. These smaller, more densely packed components, 
by virtue of the heat dissipation within a small volume, frequently operate at exces- 
sively high temperatures. These high temperatures result in component derating, 
performance degradation and accelerated failure. [Ref 1] 

Successful thermal management of electronic systems, under development in 
the latter part of the 1990s will require the removal of as much as 500 W from a 
single chip, at heat fluxes between 50 to 100 W/cm?. Volumetric heat release rates 
can also be expected to increase dramatically and are likely to exceed 10 W/cm*. 
Silicon chips, with embedded bipolar circuits, have traditionally been maintained at 
temperatures ranging from 65° C, for commercial computers, to between 110° C and 
125° C for military equipment. [Ref 2] 

Due to necessity, military electronics are frequently operated in confined 
spaces aboard ships and aircraft. This limited space impacts significantly on the 


space available for extensive cooling equipment. Air-cooled avionics systems are 


currently the most frequently encountered cooling systems used in military appli- 
cations. However, it 1s common practice with military avionics to upgrade the 
operational capability of a ship or aircraft by upgrading only one or two pieces of 
avionics equipment at a time. With an air cooled system, this is very ineffective. 
By replacing equipment in one of these systems with equipment of a different size 
and shape, the cooling airflow of the original design is disturbed and the resulting 
airflow around the new (as well as the old) equipment is less than optimal. The 
original cooling system was designed for equipment configured in a certain manner 
and little attention is paid to the reoptimization of the cooling system when changes 
are made. Moreover, cost restrictions are also an important factor in the design of 
military cooling systems. 

One solution for the military applications problem is the development of a 
universal rack for shipboard or aviation use. If this rack were modular, it would allow 
for a great amount of flexibility in design and provide significant cost savings. With 
a modular design, upgrades to equipment could also be accompanied by an upgrade 
to the cooling rack itself with very little additional cost or difhculty. By developing 
one universal rack for al] military applications, a significant improvement in cooling 
systems for updated designs can also be realized. This could be accompanied by 
additional cost savings through the elimination of unique cooling systems for every 
different avionics suite. 

The rack would be a structure that could accomodate the placement of several 
modular electronic packages. One might imagine a “shoeboz” configuration with, say, 
16 “holes” in which 16 electronic packages could be inserted. Each of the packages is 
presumed to contain electronic components mounted in a variety of ways (i.e printed 
circuit boards or on brackets attached to the walls of the package). Such a structure 


would possess at least 64 flow passages exclusive of flow passages that are part of 


the fluid source or pump. The objective is to assure that each flow passage carries 
enough fluid to absorb the dissipated heat that is somehow conducted to the rack 
structure. Of course, a heat exchanger may be required to transfer the heat to the 
ultimate heat sink which may be the environmental! air or sea water. 

Thus, the basic assumption is made that a modular, water cooled avionics 
rack can provide sufficient cooling for any piece or combination of avionics equipment 
if enough water flow paths are used, the water is at the appropriate temperature and 
the water is properly distributed within the passages. To determine if the cooling 
medium, water, is being sufficiently distributed within a modular cooling rack, an 
analysis of the flow and pressure distribution of the water is required. This can be 
accomplished using a computer code utilizing node analysis. Such a computer code, 
allowing for a variable number of branches and junctions, is presented here as a first 


step in the development of a universal cooling rack for military applications. 


C. BASIC THEORY AND APPROACH 
1. Node Analysis 

The analysis of the cooling rack is based on the stipulation that any 
size (diameter and length) of passage may be used in the construction of the rack. 
Variables in the program include the density and viscosity of the water, the number 
and location of flow sources, the number of water paths entering each junction, the 
shape of the flow passages (circular or rectangular) and the length and equivalent 
diameter of each passage. The purpose of this degree of flexibility is to allow for easy 
redesign of a rack in the event that it does not meet the requirements which are the 
proper distribution of cooling water in each section of pipe. By varying either the 
rack configuration or the state of the water via computer input, a rack that provides 


the proper flow distribution can eventually be proposed. 


The calculations for the flow in the passages employs a matrix oriented 
procedure used in network analysis. The network analysis approach can be tailored 
to flow in passages by proposing an analogy of the current” to fluid flow and the 
”voltage” to the fluid pressure. The ”resistance” is then attributed to the friction in 
each length of pipe. Therefore, each length of pipe will have a resistance associated 
with it, and possibly a pressure source as well, depending on the design of the cooling 
rack. 

2. Laminar Flow 

The computer code is designed to calculate a laminar flow distribution. 
If the flow is in transition or turbulent, there is a significant increase in the amount 
of frictional resistance. For Reynolds numbers (Re) less than or equal to 2100, the 
flow is considered laminar. For Reynolds number between 2100 and 10000, the flow 
is in transition and the flow is turbulent if the Reynolds number exceeds 10000. The 
computer code calculates the Reynolds number for each flow passage and provides 


a warning if the Reynolds number exceeds 2100. 


De sCOrr 


e Chapter II explains and details the basic code required to analyze the flow 


and pressure distribution of water in a proposed cooling rack design. 


e Chapter III describes the computer code, its essential capabilities and limita- 


tions, associated subroutines, input requirements and final output. 


e Chapter IV presents the results from several case runs that exhibit the flexi- 


bility and capability of the method. 


e Chapter V concludes with future development efforts and the application po- 


tential of this computer code. 


II. AVIONICS COOLING RACK PROBLEM 
FORMULATION 


A. GENERAL 

The flow distribution of the fluid (in all branches) of the cooling rack is 
determined using a matrix oriented solution technique. The general equations and 
the matrix solution are presented in what follows. 

1. COOLING RACK BASIC STRUCTURE 

The basic structure of the cooling rack is a series of fluid passages fit- 

ted together using standard junction components. An almost unlimited number of 
passages may be fitted together. An example of one possible simple combination of 


passages with one flow source is presented in Figure 2.1. 


@ 
low 


L 
ne, 





Figure 2.1: Example of cooling rack design. 


This example illustrates how variable sized flow passages may be used. 


The rack is composed of 6 straight passages called branches and n; junctions called 


nodes. The branches are numbered because the computer code requires that basic 
information for each flow passage branch be entered separately. The nodes are also 
numbered and those are shown in circles. Note also that there is a pressure source 
located at branch-1. For convenience, the numbering sequence in Figure 2.1, begins 
with a 1 on the upper left side and proceeds to the right and then down. Although 
this numbering sequence is arbitrary, maintaining a consistent technique facilitates 


later adaptations to the rack and comparison with other rack configurations. 


B. VARIABLES IN DESIGN 


As stated in a previous section, the length, shape and the diameters of the flow 
passages may vary. The code is written for either circular or rectangular passages. 
For rectangular passages an effective diameter is calculated. Although the actual 
diameter or effective diameter may vary within one rack design, it is assumed that 
either all circular or all rectangular passages are used for any one rack configuration. 
All length and diameter size information is entered in the initial part of the program 
as bx 1 matrices named ELL(IB) and D(IB), respectively. The density and the 
viscosity of the cooling water are input as the variables RHO and MU. The other 
input and output variables used in the computer code are summarized, for the 


reader’s convenience, in Table 1. 


C. DEVELOPMENT AND LINEARIZATION OF 
FLOW EQUATIONS 


1. Laminar Equations 
The basic equation used to determine the change in pressure or pressure 


loss within a fluid passage is derived from the D’Arcy equation 


_ f LV? 


in = 
29d, 





(2.1) 


TABLE 2.1: INPUT AND OUTPUT VARIABLES (FORTRAN) 


Variable 


Al 
B 
Bl 
D 
ELL 
IB 
MU 
N 
NF 
NT 
PS 
RHO 


Explanation 


Cross sectional height 
Number of branches 
Cross sectional width 
Branch diameter vector 
Branch length vector 
Counter for branches, IB = 1 to B 
Viscosity 

Number of nodes 

Node “from” array 
Node “to” array 
Pressure source vector 
Density 


In this equation, hy is the pressure loss due to friction and d, is the 


equivalent diameter. The D’Arcy equation is valid for steady flow within passages 


running full of liquid. In eq (2.1), f is a friction factor. If the flow is laminar, the 


friction factor, f, as shown in Figure 2.2, can be represented by the equation 


64 
J= 


The Reynolds number, by definition is 


Re= pVide 


Therefore, f, is seen to have an alternate form 


64 
f= pVd. 





Substitution of eq (2.4) into eq (2.1) yields 


_ 32LV yp 
pg? 
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(2.2) 


(2.3) 


(2.4) 


(2.5) 


Friction Factor, f 
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Figure 2.2: Moody’s chart for the friction factor. 


By definition, the flow rate, Q, is 
al (2.6) 


The area, A, of a circular flow passage is, of course, given by 


d? 
A = saere (2.7a) 


but for rectangular passages of side dimensions a and 8, it is 

A— ap (2.7b) 
and for square passages where a = 8, it is 

A=a PAKS) 


For circular passages, d., is simply the passage diameter, d, and for rectangular 


passages, d, is defined as 


de 


II 


CA 

le 

where A is the passage flow area and P, for channels flowing full, is the passage 
wetted perimeter. In order to make the equivalent diameter for a circular passage 
equal to the actual diameter, d,C = 4 so that for a rectangular passage with side 


dimensions, a and 8, the equivalent diameter becomes 


2ab 
a+b 





d, = (2.8a) 


In the event that the passage is square (a = 5), the equivalent diameter becomes 


d. =a (2.8b) 


Substitution of eqs (2.7) and (2.8) into equation (2.5) yields for the 


circular passage 
ie 128LQ pu 


Te (2.9a) 
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for the rectangular passage 


_ 8LuQ(a + b)? 








a 2.96 
‘~~ pg(ab) = 
and for the square passage 
32 ae. 
pe 
f re (2.9c) 
To simplify eq (2.9a) further, let 
oe 
and 
ei 
i aie (2.10) 
This yields the equation 
ie Hl) 
which is of the same form as 
Yel 


Similar adjustments can be made to eqs (2.9b) and (2.9c). When these adjustments 


are made, it is seen that for a rectangular passage, eq (2.9b) gives 
_ SLy(a + 6)? 
~— pg(ab)> 
and that for a square passage eq (2.9c) provides 
R= ae 
pga 


Additional losses (other than passage friction already discussed) occur 


R (2.11) 





in flow passage systems and cannot be disregarded without appreciable error. Com- 
pensation for entrance and exit losses must be considered in the cooling rack system 
by adding an equivalent length of straight pipe. The equivalent length used for 
square edged entries is 20 diameters (or equivalent diameters) and for exits, the 
equivalent length is 40 diameters (or equivalent diameters) [Ref. 2]. Therefore, the 


additional length of 60 passage diameters accounts for both exit and entrance losses. 


10 





D. NUMERICAL SOLUTION SCHEME 

We first consider a general branch (the k"* branch) as shown in Figure 2.3. 
The branch can contain a pump with pressure head, p,,, with an associated section 
of pipe of resistance, Ry. The branch carries a flow, g,, and possesses a total head 


loss, pz. Notice that the orientation of the head loss is in the direction of the flow 


(from “+” to “—”). 





Figure 2.3: Branch with pressure source. 


Also observe the orientation of the flow, g,, and notice that by the analogy 
to Ohm’s law and because of compatibility (the sum of the losses must match the 


head available) 


Pr = PRaqe + Dok 


This may be written for all k branches in matrix form 
P=RQ+P, 


Consider a pump and a length of discharge piping. The pump can be repre- 


sented as a flow source or as a pressure source as indicated in Figure 2.4. 


1] 





Figure 2.4: Alternative source arrangements for the development of the 
Flow source «+ Pressure source transformation. 


In these figures, the R’s represent the resistance to flow of the discharge 
piping and, in the case of the pressure source, the + and — indicate the discharge 


(+) and suction (—). In Figure 2.4a, continuity dictates that, 


gz tq 
Hs i 
or 
p= q+ Aas 


In Figure 2.4b, compatability (the pressure drops around the simple fluid loop must 


match) 
p= iyps 


If there is to be an equivalence between the two (Figures 2.4a and 2.4b) then 
R=R,=R, 


and 


= Rqs 


i 


or indeed 
q,= 
ae 


Thus, when we represent a pump as a pressure source with a resistance, we can 
immediately transform this to an equivalent flow source. The computer code 1s 
written for the source input to be a flow source. 


Note also that continuity dictates that for the k’ branch 
dk = YePk + Qsk 
And for all k branches this can be written in matrix form as 
Q=YP+Q, 


Next consider a rack containing 6 branches and n,; nodes. The n‘* node is the 
datum node (the node at the suction end of the pump), and we may set down an 
n, x 6 augmented node-branch incidence matrix A® with elements 


+1 if branch 7 leaves node z 
—1 if branch 7 enters node z 
0 if branch 7 is not incident 
upon node 2 


ay = (2.10) 
For example in the network displayed in Figure 2.5 with its oriented graph shown 
in Figure 2.6, there are n, = 5 nodes and 6 = 6 branches. For this network, the 


augmented node branch incidence matrix A® is 


eo tO Ue Sai 

0 -l 1 Sa 0°70 

A° = 0 O -l Le tc 
ees) eae et 

ale SOO 0 20" =0 


Note that for each column, which represents a single branch, there should 


only be two non-zero entries, a 1 and a —1. These entries correspond to the nodes 
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Figure 2.5: Network for example. 





Figure 2.6: Oriented graph for network in Figure 2.5. 


at the extremities of the branch. A branch that leaves a node shows a +1 while a 
branch that leaves a node shows a —1. 

Next, we form a node-branch incidence matrix which is n x 6. This is done 
by deleting the row in A® that corresponds to the datum node. In our example, 


row-5 is the datum node and thus 


leet) Same (0 () 
0-1 1 #21 = #O 90 
ae Ve len US ol 0 
Pe 


The vector Q is defined as a 6X1 vector representing the flow in all & branches 
and its elements are designated as q,. The product of A with Q, AQ, represents a 
series of equations, one equation for each node. These equations represent the sum 
of all flow in and out of the respective node. Continuity demands that product, AQ 
be null 
AQ=0 
and this simple matrix equation is a statement of continuity for the whole system. 


To illustrate 


q1 
l J UES) q2 (qi + q) 0 
DS ee Oa ae (—q3 + qs) 0 
SO i US Ge (—9s — Gs + 96) 0 
d6 


The pressures or heads at the nodes are represented by an n x 1 vector 
designated by H and the branch pressure drops are represented by an n x 1 vector 


designated by P. Reference to Figure 2.6 shows that 


P2 = hi-ho 
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P3 = ha—hz 


Pa = hoa—hg 
es — h3 — hg 
lo hy 


1 0 0 0 
1 —-l 0 0 hy 
= _@ig@d 1 —-1 0 h» 
b= aE 0 ] 0 —-1 hg 
0 0 1 —-1 hs 

| 0 0 0 Ty 
and it can be observed that 
C=A? 


Refer now to the general branch in Figure 2.3 and the vector representation 
for the flow in all & branches. The 6 x 1 vector, Q which represents the branch flow 


rates has been shown to be 


Oi =MP-4e. (2.11) 


Again, we carefully note the orientation of the flow source. Here Y is an n x 


nm source admittance matrix. Premultiply eq (2.11) by A to obtain 
AQ=AYP+AQ, 

But, P = A?H and AQ = 0, thus 
AYA'TH + AQ, = 0 


The node equations are then formulated 


where 


Y, =AYA! 


and 


Q = -AQ, 


Therefore, the node pressures H, can now be determine by 
H=Y,'Q 


Although many techniques exist for the solution of large sets of linear si- 
multaneous algebraic equations, the most efficient computationally, appears to be 
the Cholesky reduction followed by back subsitution that is employed in the Gauss 
elimination method. The only restriction on the use of Cholesky reduction is that 
its use is confined to symmetric, positive definite matrices. Here it is fortunate that 
Yn is always positive definite. 

The Cholesky reduction is based on the premise that there is a unique lower 
triangular matrix that permits a factorization in the form of A = LL? if the matrix 
A is symmetric and positive definite. Computationally, the Cholesky reduction is a 
very efficient technique in that it only requires n(n + 1)/2 storage locations for L, 
rather than the usual n* locations required by other methods. [Ref. 3]. The number 
of operations using the Cholesky reduction is approximately n°/6 rather than the 
usual n°/3 required for most other decompositions (Ref. 3]. 


The Cholesky method to solve the basic system 
AX=B 
is based on finding the solution to two equivalent systems: 
L'C=B and LX=C 
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The elements of C are determined by the algorithms 


b, 
qQ = 
$11 


and 


1—] 
b —S— sace 


S43 


Once C is known, X can be found using back substitution as employed in the Gauss 


elimination method, [Ref. 5] that is 





and 


III. DESCRIPTION OF THE COMPUTER 
CODE 


A. PROGRAM STRUCTURES AND CAPABILITIES 


1. Restrictions and Limitations 
The computer system used is an IBM or IBM compatible personal com- 
puter. The current program fixes the maximum number of individual fluid passages 
allowed at 300. For proper use of this program, the computer must have a minimum 
storage requirement of 2 M Bytes. A detailed computing time study for the program 
has not been undertaken because the computing time changes exponentially with 
the number of nodes and branches selected in the proposed cooling rack design. 


The computer code can be operated using the following information 


e Laminar Flow Conditions 
e Density of water varying from x to y 


e Viscosity of water varying from x to y 


Up to 100 branches 

e Up to 100 nodes 

e Up to 40 sources 

e Any diameter (or equivalent diameter) flow passages 


e Branches of any length 


Metric or English units 


1g 


2. Structure of the Main Computer Code 
The the main computer code including the three associated subroutines 
are in Appendix A. The main code is essentially divided into two major sections, 
initialization/input and laminar flow solution with its associated output. The main 


function of the initialization/input section are: 


1. Set up the problem formulation by reading in the necessary values (node from, 
node to, length and cross sectional dimensions) for each branch of the cooling 


rack. 
2. Read in the viscosity and density values. 


3. Read in location and characteristics of each pressure source. 


The second section calculates the pressure and flow distribution in the 
rack. It includes the subroutines MATMUL, DECOMP and CHOLESKY. The main 


functions of this section are to 


1. Develop the node-branch incidence matrix, A, using the information read in 


section one. 
2. Calculate the effective branch lengths to account for losses at the nodes. 


3. Calculate the area and effective diameter of each flow passage cross sectional 


area. 
4. Calculate the branch admittance matrix, Y. 
5. Develop the flow source matrix, Qs, from information read in section one. 


6. Calculate the transpose of the matrix A. 
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10. 


Ne 


12. 


13. 


14, 


B. 


. Calculate the node admittance matrix, Yn. 
. Calculate the node flow source vector, Is. 


. Calculate the node pressure vector, H. 


Calculate the inverse of Yn. 

Calculate the branch pressure vector, P. 
Calculate the branch flow rate vector, Q. 
Calculate the branch Reynolds numbers. 


Provide readouts of all branch flow rates, all node pressures and all branch 


Reynolds numbers. 


DESCRIPTION OF SUBROUTINES 
1. Subroutine MATMUL 


This subroutine (Ref. 5] multiplies an n x m matrix by any m x £ matrix 


to form an n x | matrix and is called whenever matrix multiplication is required. 


2. Subroutine DECOMP 


This subroutine [Ref. 6] performs the decomposition of the symmetric, 


positive definite matrix, Y, using the Cholesky reduction. It is called by subroutine 


CHOLESKY. 


3. Subroutine CHOLESKY 


This subroutine [Ref. 7] provides the solution of a linear system of 


equations using the Cholesky decomposition method for positive definite matrices. 


This subroutine calls subroutine DECOMP. 


Zl 


IV. RESULTS AND DISCUSSION OF CASE 
RUNS 


A. CASE RUN ONE 


Case one was run using the configuration shown in Figure 2.1 and the follow- 


ing data: 


1. Water density = 62.4 lb/ft* 


bo 


. Water viscosity = 8.8 x 1074 lb/sec-ft 

3. Basic shoebox design with 8 nodes and 12 branches 
4. Length of end branches (1-4 and 9-12) = 2.2 ft 

5. Length of middle branches (5-8) = 3.0 ft 

6. Circular passages with a diameter of 0.3 inches 

7. One flow source at branch one 


8. Strength of source was 0.4 gal/min 


The results from case one are in Appendix B. It should be noted that the 
Reynolds number in each branch indicates laminar flow and that there is cooling 
water in all sections of the sample rack design. These results demonstrate that 
the analysis technique presented here, does allow for the determination of whether 
or not a proposed cooling rack design will indeed provide a proper distribution of 
cooling water. It also verifies that the laminar flow assumption is valid for some 


design parameters. 


tO 
tO 


B. CASE RUN TWO 
Case run two demonstrates the flexibility of the program through a varia- 
tion in some of the input parameters. The basic configuration is still the same 


configuration shown in Figure 2.1. 


1. Water density = 999.5 kg/m*° 

2. Water viscosity = 13.1 x 1074 kg/m-sec 

3. Basic shoebox design with 8 nodes and 12 branches 

4. Length of end branches (1-4 and 9-12) = 1.0m 

5. Length of middle branches (5-8) = 1.2 m 

6. Rectangular passages with a height of 1.3 cm and width of 1.4 cm. 
7. One flow source at branch one 


8. Strength of source was 1.8 lit/min 


The results from case two are located in Appendix B. 


V. CONCLUSION 
A. GENERAL COMMENTS 


This computer code has been developed as an initial step in the total design 
of a modular cooling rack for avionics equipment. In itself, the code details a 
specific design technique and allows for the determination of whether one proposed 
configuration, including source location, characteristics of the cooling water, and the 
size and shape of the proposed flow passages will indeed provide a proper distribution 
of the cooling water. Without proper coolant distribution, the cooling rack will be 


inefficient and perhaps, totally ineffective. 


B. ENHANCING THE MAIN PROGRAM 
CAPABILITY 


The next step in the development of a computer code to provide more flex- 
ibility and range in the rack design is the inclusion of the turbulent flow regime 
within the coolant passages. After turbulent flow is effectively incorporated, the 
final phase in the total rack design is a modification for heat input in the individual 


coolant flow passages. 
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APPENDIX A 


PROGRAM COOL 
Cooling Rack Calculations Using Node Analysis 
INITIALIZATION 


MU: fluid viscosity (lb/ft-sec) 
RHO: fluid density (lb/f£t*3) 
N: number of junctions in system 
Ni: number of junction minus the datum node 
B: number of branches in system 
G: acceleration of gravity (ft/sec”2) 
PI: pi 


Cl] Ub sd AOS 42) ed 2) ed) eae 


TNTVEGE Rede i So lj, hk, Ob cS 

INTEGER NT(40) ,NF(40) ,X, IUNITS 

REAL MU,RHO,G,PI,A(40,40) ,D(40) ,ELL(40) ,R(40) ,MU1 
REAL QS(40,1) ,AT(40,40) ,QS1,11(40,1) , LENGTH 

REAL E(40,1),P(40,1),Y(40,40) , YNI (40,40) 

REAL YP(40,1),Q(40,1),AY(40,40) ,YN(40, 40) 

REAL V(40,1),RE(40,1),a1(40) ,b1 (40) ,AREA(40) 

REA IS(407 1) sit 40,1) ,RHOL 

CHARACTER ANS, ANS1,ANS2 


G=32.174 
PI=3 21415926 
IUNITS=0 


KKK KKK KKK KKK KKK KKK KKK KKK KK KK KKK KKK KKK KKK Ka KKK KKK KKK KKK KK KKK KKK KKK Keka Kk 


THIS PROGRAM DETERMINES THE FLOW RATE AND PRESSURE DISTRIBUTION 
OF COOLING WATER WITHIN A VARIABLE-SIZED AVIONICS COOLING RACK. 
VARIABLES INCLUDE RACK DIMENSIONS AND WATER CHARACTERISTICS. 


KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KEK KKK KKK KEKKKKaKa KKK KK KKK kKaKk 


<p (ages Ge EG A GS oo anes Va We) Pie) (et Go Jew 


1040 FORMAT (///' ENTER THE TOTAL NUMBER OF NODES IN THE RACK SYSTEM 
tt ey oN) 

1042 FORMAT (BN, I13) 

1043 FORMAT (//' ENTER THE TOTAL NUMBER OF BRANCHES IN THE RACK SYSTEM 
Be Palio aN 


a | 
— 


bee WRITE (IOT, 1040) 
READ (IN,1042) N 
WRITE(IOT,1045) N 
IF(N .LT. 0 .OR. N .GT, 200) eh 
WRITE (IOT, 6565) 
GOBYO 222 
BN Dae 
PLO WRITE (TOT; 1043) 
READIN, 1042 )3 
WRITE(IOT,1046) B 
IF(B .LT. 0 .OR. B .GT. 100 wrt 
26 


WREEE (TOT, 6567) 


CO nO 223 
END IF 
045 FORMAT (//' YOU ENTERED THE NUMBER OF NODES AS:’ ,1X,13\) 
046 FORMAT (//' YOU ENTERED THE NUMBER OF BRANCHES AS:’ ,1X,13,/,\) 
565 FORMAT (/' THE NUMBER OF NODES MUST BE GREATER THAN ZERO AND LESS’, 
Pax, THAN 100.7, \) 
567 FORMAT (/' THE NUMBER OF BRANCHES MUST BE GREATER THAN ZERO AND LESS 
+ THAN 100.’,\) 
N1=N-1 
DO 5 I=1,B 
OS(I,1)=0. 
5 CONTINUE 
5 CONTINUE 


Welter (TOT, 7676) 
576 FORMAT (/’ Do you want to work in British units or SI units (B/S)? 
= NY 
READ (IN, 7677) ANS 
5/7] FORMAT (A1) 
IF(ANS .EQ. ‘B’ .OR. ANS .EQ. ‘b’) GO TO 7800 
MeeeNS .EO. ‘S’ .OR. ANS .EQ. ’s’) GO TO 7678 
GO TO 55 
300 IUNITS = 1 


WRITE (IOT, 1000) 

900 FORMAT (//’ Input Viscosity ( x 10°4 lbm/ft-sec)',2xX,\) 
READ (IN,1003) MU 
MU1=MU* .0001 
WRITE (IOT,1002) 

902 FORMAT (/’ Input density (lb/ft%*3)’,2x,\) 
READ (IN,1001) RHO 

901 FORMAT (BN, F8.4) 

903 FORMAT (BN, F7.6) 

904 FORMAT (BN, 13) 
Ge TO 7801 


578 WRITE (IOT, 7602) 
502 FORMAT (//' Input Viscosity (Kg/m-sec x 10°4)’,2X,\) 
READ(IN,1003) MU 
WRITE (IOT, 7603) 
503 FORMAT (/’ Input Density (Kg/m*3)’,2X,\) 
READ (IN,1001) RHO1 
MU1=MU* .00006719 
RHO=RHO1* .06243 


A: NODE-BRANCH INCIDENCE MATRIX 
a7) 


C INITIALIZE A MATRIX 
Cc 
7801 DO 20 I=1,N1 

DOr2 07 U=_ 75 


20 A(I,J)=0. 
C 
DO 22 talc 
NF (I) =0 

22 NT (I) =0 

Cc 

C 
‘e DETERMININATION OF PASSAGE SHAPE 
C 
C 
Cc 

8000 WRITE (IOT, 8001) 

8001 FORMAT (/' IF FLOW PASSAGES ARE CIRCULAR, ENTER A ZERO.’ ,/, 

+ ! IF PASSAGES ARE RECTANGULAR, ENTER A ONE.’, 2X, \) 


READ(IN, 1042) X 
[E(x EO. 0) HEN 
GO TO 8004 
END IF 
TE (XxX 2EOy 2) aaneEN 
GO TO 8005 
END IF 
8003 WRITE (IOT, 8006) 
8006 FORMAT (/' ERRONEOUS INPUT’ \) 
GO TO 8000 


INITIAL DATA INPUT TO DEVELOP A,ELL AND D MATRICES FOR 
RECTANGULAR PASSAGES 


ELL: LENGTH 
D: DIAMETER 


9QaqaaQagaiaQaQqggaaa 


8005 CALL CLS 
IF (IUNITS .EQ. 0) THEN 
DO 1493 J=1,B 
WRITE(IOT,1051) J 
WRITE (IOT, 1492) 
1492 FORMAT (/’ FROM NODE, TO NODE, LENGTH(m), HEIGHT(cm), WIDTH(cm)’, 
+ /) 
READ(IN, 8202) NF(J),NT(J),ELL(J) ,al(J) ,b1 (J) 
1493 CONTINUE 
7272 FORMAT (/’ YOU ENTERED THE FOLLOWING DATA.’ \) 
7273 FORMAT (/’ BRANCH FROM NODE TO NODE LENGTH(m) HEIGHT (cm) ’, 2X, 
+ “WiDUH (em jee) 
7274 FORMAT (2X,13,6X,13,8X,13,2X,F9.5,4X,F7.5,4X,F7.5) 
7275 FORMAT(/’ IS THIS CORRECT (Y/N)? ‘\) 
CA, Gls 
16 WRITE (IOT, 7272) 
WRITE (IOT, 7273) 
DOV7277 J=l17s 
WRITE (IOT, 7274) J,NF(J) ,NT (J) ,ELE) aii Pei 
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e/'/] CONTINUE 
276 Weer (1OT,7275) 
READ (IN, 7677) ANS 
Maas .hO., ’Y’ .OR. ANS :EQ.’y’) GO TO 7278 
mens 2-6. ’N’ .OR. ANS .EQO. ‘n’) GO TO 7280 
eeero 7276 
278 ee 7279 J=1,B 
Pm) =Bob(d) *3.28 
eee) =al(J)*.3937 
fem) =b1 (J) *.3937 
B79 CONTINUE 
GO TO 1495 
280 cong, I NUE 
oo 1 FORMAT (/’ INPUT ONE BRANCH NUMBER TO CHANGE VALUES.’ \) 
284 Witte (IOT, 7281) 
READ (IN,1042) Kl 
282 FORMAT (/’' INPUT VALUES FOR BRANCH:’ ,1X,13,\) 
WeITe(IOT,7282) Kl 
WRITE (IOT, 1492) 
eee IN, 8202) NF (K1),NT(K1) ,ELL(K1) ,al(K1),b1(K1) 
283 FORMAT(/’ DO YOU HAVE ANYMORE CHANGES (Y/N) ?’ \) 
a > WEeene {1OT, 7283) 
READ (IN, 7677) ANS1 
Meo! .FO. ‘’Y’ .OR. ANS] .EQ. ‘y’) GO TO 16 
jeeNol) .BO. ’N’ .OR. ANS1 .EQ. ‘n’) GO TO 7278 
@2)TO 7285 
END IF 


CONTINUE 


DO 8200 I=1,B 
WRITE (IOT,1051) I 
WRITE (IOT, 8201) 
101 FORMAT(/’ FROM NODE, TO NODE, LENGTH(FT), HEIGHT(IN), WIDTH (IN) 
+’,/) 
102 FORMAT (213,F7.4,2F6.4) 
MeamiiINn,8202) NF(I),NT(1I),ELL(I),al(1I),b1(I1) 
100 CONTINUE 


eon, CLS 
; WRITE (I1OT, 7272) 
WRITE (IOT, 8207) 
107 FORMAT (/’ BRANCH FROM NODE TO NODE LENGTH(FT) MHEIGHT(IN)’,2X, 
WIDTH (IN) ’,/) 
174 FORMAT (2X, 13,7X,13,7X,13,4X,F9.5,6X,F7.5,3X,F7.5) 
DO 8277 J=1,B 
Wie TOT, 7374) J,NF(J) ,NT (J) ,ELL(J),al (a) ,b1 (J) 
\77 CONTINUE 
'76 WRITE (IOT, 7275) 
READ (IN, 7677) ANS 
IF(ANS .EQ. ‘’Y’ .OR. ANS .EQ.’y’) GO TO 1495 
TF(ANS .EO. ‘'N’ .OR. ANS .EQ. ‘n’) GO TO 8280 
| GO TO 8276 
80 CONTINUE 
84 WRITE (IOT, 7281) 


+ 
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READ (IN,1042) K1 

Wiekos (lOT 7262) een 

WRITE (1Om 72 oir) 

READ (IN, 8202) NF (K1) ,NT(K1) , ELL(K1) /ait ke ei 
B2Zzie5 WRITE (OT) 7 2cc0 

READ (IN,7677) ANS1 

IF (ANSI .EQ. ‘Y’ .OR. ANS] .EQ. “y" )Gomuomens 

IF(ANSI .EQ. ‘'N’ .OR. ANS1 .EO. ’n’)) GOu nema 

GOTTOReZeS 
HAS DO 1496 I=1,B 

LE CE pe Ny) ren 


re NES Gere) eal 
END IF 

IF (NT(I) .LT. N) THEN 
PANG ck) = 7 i 
END IF 


ail(I)=al1(I)/12 
Hie oil) 7/12 
AREA (I) =al (I) *b1 (I) 
D(I) =(2* (al (1) *b1 (1) )) / (al sbi 
LENGTH=ELL (I) 
ADDL=60*d (I) 
ELL(I)=LENGTH + ADDL 

1496 CONTINUE 

GO TO 8300 


RECEIVE DATA FROM KEYBOARD TO DEVELOP A,D AND ELL MATRICES FOR 
CIRCULAR PASSAGES 


I 6ey9) FORMAT (1X,13,8X,13,5X,F7.4,4X, F6.4) 
2050 FORMAT (213, 2F7.4) 
8004 CALL CLS 

TP (UNS. -FO. 0), THEN 


COETC3s 7. 
ELSE 
DO 2 ee ee 
OS a FORMAT (/’ INPUT THE FOLLOWING DATA, SEPARATED BY COMMAS FOR, 1X, 
+ BRANCH:’,13,2X, \) 
HOS 3 FORMAT (/’ FROM NODE,TO NODE, LENGTH (FT) , DIAMETER(IN) ', 2X, /) 


WR en On LO S12 1 
Whee (Om Le S33 
READ (IN, 9050) NF(I),NT(1), ELL( hye 


25 CONTINUE 
Gn, CEs 
18 Weiler hom, 7oa2) 
WRITE (IOT,1059) 
1059 FORMAT(/’ BRANCH FROM NODE TO NODE LENGTH(FT) DIAMETER(IN)’, 
+ /) 


DOS 2 77 ao— i> 

WRITE (1OT, 1155) J,NF(d) NT (9) (ELL Gea 
aA Sys FORMAT (2X,13,9X,13,8X,13,6X,F/.4,4x% 86 =) 
Zo 5 FORMAT (2X,13,7X,13,7X%,13,5X,F 7 473 59s ee 
S207 CONTINUE 
S276 WRITE (1OT, 7275) 

READ (IN,7677) ANS 
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PeaicometOrm y mmOR. ANS -EO.’y’) GO TO 26 
PeenNcmemene  N §.OR. ANS .EO. ‘n’)™GO° TO 9280 
GOTO 9276 

280 CONTINUE 

284 Ricehe (lor ,7281) 
READ (IN,1042) Kl 
Were (1OT,7282) Ki 
Wie (TOT, 1053) 
mee (IN, 9050) NF(K1) ,NT(K1),ELL(K1) ,D(K1) 

285 Weeris (LOT, 7283) 
READ (IN, 7677) ANS1 
jeeaniot .BO. ‘Y’ .OR. ANS1 .EQ. ‘y’) GO TO 18 
[iewaNS] .BO. ’N’ .OR. ANS1 .EQ. ‘’n’) GO TO 26 
cee 9285 
END IF 
Ge TO 26 


73 ber29 J=1,B 
MiereE(lOT,1051) J 

5 FORMAT (/’ FROM NODE, TO NODE, LENGTH (M) , DIAMETER (CM) ’, /) 
ete (TOT, 1055) 
meen IN,9050) NF(J) ,NT(J) ,ELL(J) ,D(J) 


] CONTINUE 
CALL CLS 
] PEs (1OT, 7272) 


rors (TOT, 1077) 
7 FORMAT (/’ BRANCH FROM NODE TO NODE LENGTH(M) DIAMETER(CM)’, 
+ /) 
DO 1277 J=1,B 
WRITE (IOT,1255) J,NF(J) ,NT(J),ELL(J) ,D(J) 
277 CONTINUE 
276 Wrere {1TOT, 7275) 
READ (IN, 7677) ANS 
MeeeNs .FO. ’Y’ .OR. ANS .EQ.‘y’) GO TO 111 
MWEEeANS .EQ. ‘’N’ .OR. ANS .EO. ‘n’) GO TO 1280 
GerTO 1276 
80 CONTINUE 
184 WRITE (IOT, 7281) 
READ(IN,1042) Kl 
WRITE (IOT, 7282) Kl 
WRITE (IOT,1055) 
READ (IN, 9050) NF(K1) ,NT(K1) , ELL (K1) , D(K1) 
385 WRITE (IOT, 7283) 
READ (IN,7677) ANS1 
IF(ANS1 .EQ. ‘Y’ .OR. ANS1 .EQ. ‘y’) GO TO 19 
/ManNst .FO. ’N’ .OR. ANS1 .EQO. ‘n’) GO TO 111 
GO TO 1285 


.1 Be 23 J=1,8 

| Bei) =ELL(J) *3.28 
Beg y=D(J)*.3937 
CONTINUE 


- »- ART. «~,* 


wers3s I=1,85 
Ba 


gaa0a0Q 


Q 


NQAaQAaQqaQqNagAnNaaANn 


33 


8300 


30 


35 


40 


OS; 


OB: 
Osi: 


dO INS 
GAG 
Oy 
2 Oey, 


IF (NF(I) .LT. N) THEN 
PSCONIS IL) ploy Weal 


END IF 

IF (NT(I) .LT. N) THEN 
PN) | Tr) == 
END IF 


Dil) =pigh ee 

AREA (I) =PI* (D(I) **2) /4 

LENGTH=ELL (I) 

ADDL=60*D (I) 

ELL (I) =LENGTH+ADDL 
CONTINUE 


DO 30 K=1,B 
R (K) = (RHO*G*AREA (K) * ((D(K) ) **2)) /(32.*MUL*ELL (XK) ) 
CONTINUE 


DO-35 J=17B 

DO 35 K=1,B 
Y¥(J,K) =0 

CONTINUE 


DO 40 I=1,B 
Nel) — ht 1) 
CONTINUE 


MATRIX OF FLOW SOURCES 


NUMBER OF FLOW SOURCES 
SOURCE BRANCH 
SOURCE STRENGTH 


FORMAT (/’ Input the number of flow sources:’ ,3X,\) 
FORMAT (//' Input the branch of source’ 73 22 
FORMAT (//' Input the strength of the source (gal/min) | 33908 
FORMAT (//’ Input the strength of the source (liter/min) :’ ,3X,/) 
WRIEE C1lOL 7 boas; 
READ(IN,1004) S 
DO w50. = s 
WRITE (IOT,1016) I 
READ (IN,1004) QB 
IF(IUNITS .EQ. 1) THEN 
WRLEE GlOl, 107) 
READ(IN,1001) QS1 
ELSE 
WRITE (IOT, 9017) 
READ(IN,1001) QS1 
OSi=O0S1"*,. 2642 
END IF 
0S1=0S1/444 
DO 60 K-i-e 
oe 


bee (kK HO. OB) THEN 
Senko) -OS 1 


ELSE 
Os(K,1)=0. 
END IF 

5 0 CONTINUE 

> 0 CONTINUE 


feeeee (TOT, 3012) (OS(1I,1),1=1,B) 


MATRIX MANIPULATION TO DETERMINE INDIVIDUAL BRANCH FLOW AND 
PRESSURE DISTRIBUTIONS 


AT: TRANSPOSE OF MATRIX A 
AY: PRODUCT OF MATRIX A AND MATRIX Y 
YN: PRODUCT OF MATRIX AY AND MATRIX AT (YN=AYAT) 


be 80 I=1,N1 
DO 80 J=1,B 
30 AT(J,1I) =A(I,J) 


CALL MATMUL(AY,A,Y,N1,B,B) 


CALL MATMUL(YN,AY,AT,N1,B,N1) 


TS NODE FLOW SOURCE VECTOR IS= -AQS 


CALL MATMUL(IS,A,QS,N1,B,1) 


B50 I=1,N1 
roel 2 )—--1S(1, 1) 
ore) = 1S (i, 1) 

10 CONTINUE 


YNI: MATRIX EQUAL TO MATRIX YN BEFORE CHOLESKI INVERSION. AFTER 
DECOMPOSITION IT HOLDS THE UPPER TRIANGULAR MATRIX. 

ISI: MATRIX EQUAL TO MATRIX IS BEFORE INVERSION. AFTER INVERSION, 
IT HOLDS THE SOLUTION TO E=YN(-1)1S 


| DO 95 I=1,N1 
DO 95 J=1,N1 
YNI (I,J) =YN(I,d) 
f CONTINUE 


53 


Noon. e200 


Q 


qgqaqaqn 


GALL CHOLESKY (YNI,IS1,N1i,N®) 


NODE VOLTAGE MATRIX E=YN(-1)IS 
BRANCH PRESSURE P=ATE 

BRANCH FLOW RATE Q=YP+QS 

PRODUCT OF MATRIX Y AND MATRIX P (YP) 


tJ 10 ‘0 BI 


DO 96 I=1,N1 
E (12S it 1) 
96 CONTINUE 


CALL MATMUL (P,AT,E,B,N1,1) 


CALL MATMUL(YP,Y,P,B,B,1) 


DO 120 I=1,B 
ONG eer OS 74) 
1210 CONTINUE 
CALL CLS 


REYNOLDS NUMBER CALCULATIONS 


DO 5000 I=1,B 
V(I,1)=ABS(Q(I, Bete ca td )) 
RE (ie) =1o* Vv (1,1) *a rim 
IF(RE(I,1) .GT. 2100) eke 
WRITE (IOT,SO0O1) I 
WRITE (IOT, 5002) 
5001 FORMAT (/’ REYNOLDS NUMBER IN BRANCH’ ,13,2x,’EXCEEDS’ , \) 
5002 FORMAT (/' 2100. THE FLOW IS NOT LAMINAR.’ 27055 
ELSE 
END IF 
5000 CONTINUE 


L500 FORMAT (3X,13,7%, F11l.8,2%, 2F1525,1 3) 0 oF 
nO Sees ae 13,5X%,F11L.8,3X%,2F15.5, 23 oer 
1501 FORMAT (///' Patel 4X,' Flow Rate (ft/sec)’,4x,'’P in "pamm 
+ Pout 2 (psi. Se ayYAD Tse 
1601 FORMAT (///* Jester cing Rate (m/sec)',5x,’P in (N/m 2) eee 
+: ? DOU ON / Me 2) Goapele osc oe) 
F(IUNITS .EQ. 1) THEN 
WRITE (IOT,1501) 
DO 7000, 1=i775 
WRITE (IOT,1500)1,0(1I,1),E((NF(1)) 1), BUNT OG 2 oe 
7000 CONTINUE 
99 FORMAT (/' WOULD YOU LIKE A PRINTOUT OF THIS TABLE (Y/N) ?’\) 
98 WRITE (IOT, 99) 
READ (IN, 7677) ANS2 
IF(ANS2 .EO. ’Y’ .OR. ANS2 .EQ. ‘y’ ) ,G@uteeey 
IF (ANS2 .EQ. ‘’N’ .OR. ANS2 .EO: ‘'n’) GOmieese 
GO TO 98 
34 


WRITE (IPR,1501) 
DO 93 I=1,B 
Paice 500)0,0(1,1),E((NF(I)),12),E((NT(I)),1),RE(I,1) 
CONTINUE 
CONTINUE 
ELSE 
DO 7470 J=1,N1 
E(J,1)=E(J,1) *47. 88 
70 CONTINUE 
WRITE (IOT,1601) 
DO 7001 I=1,B 
Q(I,1)=QO(I,1)*.3048 
Peee(TOT,1500) 1,0(1I,1),E((NF(I)),1),E((NT(I)),1),RE(I,1) 
O01 CONTINUE 
| WRITE (IOT, 99) 
READ (IN, 7677) ANS2 
ME ANS2 .EFQO. “¥’ .OR. ANS2 .EQ. ‘y’) GO TO 47 
TF(ANS2 .EQO. ‘'N’ .OR. ANS2 .EQ. ‘n’) GO TO 82 
GO TO 48 
: WRITE (IPR,1601) 
DO 43 I=1,B 
Pim (TPR,1500)1,0(1I,1),E((NF(I)),1),E( (NT(1I)),1),RE(I,1) 
2 CONTINUE 
: CONTINUE 
END IF 
END 
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SUBROUTINE CHOLESKY (A,B,N,NX) 


SOLUTION OF LINEAR SYSTEMS OF EQUATIONS BY THE CHOLESKI 
METHOD FOR SYMMETRIC POSITIVE DEFINITIVE MATRICES. 


A: 
BB: 
Ge 
N: 
NX: 


ARRAY CONTAINING THE SYSTEM MATRIX (AX=B) 

ARRAY CONTAINING THE VECTOR IF INDEPENDENT COEFFICIENTS 
AUXILIARY VECTOR 

ORDER OF A 

ROW AND COLUMN DIMENSION OF A 


DOUBLE PRECISION B(NX,1) 
DIMENSION A(40,40) 


COMPUTE UPPER TRIANGULAR MATRIX FROM A AND STORE ALSO IN A 


CALL DECOMP (A,N, NX) 


COMPUTE THE C VECTOR AND STORE IN ARRAY B 


LOO 


1S Or 


102 


5 
20 


Bl 1B ee) 
WRITE (IOT,100) A(1,1) 
FORMAT (F18.15) 
DO 10 I=2,N 
WReeEE (Om. 1 O04) bN 
FORMAT (213) 
Write LOU eLO2 ) At igi A Chien) 
FORMAT (2F13.11) 


D=B(I,1) 

I1l=I-1 

DOr S: =a 
D=D-A(L,I1)*B(L,1) 
Bi 71) =87 e011) 
B(N,1)=B(N,1) /A(N,N) 


c COMPUTE THE SYSTEM UNKNOWNS AND STORE IN ARRAY B 


C 


i> 


2,0 
30 


N1i=N-1 
DO: 30. b=17Na 
K=N-L 
WRITE (IOT,105) A(K,K),K 
FORMAT (F13.11,13) 
K1l=K+1 


BO 2070 =KieN 

B(K oo) )=Biiea) At KJ) =3 (jy) 
RETURN 

END 


SUBROUTINE MATMUL(C,A,B,N,M,L) 
This subroutine computes the matrix operation C =A * B 


N: NUMBER OF ROWS IN MATRIX A AND C 
M: NUMBER OF COLUMNS IN A AND ROWS IN B 
L: NUMBER OF COLUMNS IN B AND C 


DIMENSION A(40,40) ,B(40,40) ,C (40,40) 
BO@ez0 T=1,N 
DO 20 J=1,L 

C(I,J)=0.0 

DO 20 K= 
,J)=C(I,J)+A(1I,K) *B(K,J) 
RETURN 
END 


eu) 


Q 


0° 0O-O°0 O70, G46 


SUBROUTINE DECOMP (A,N,NX) 


THIS SUBROUTINE PERFORMS THE DECOMPOSITION OF AS eeeittt = 
DEFINITE, SYMMETRIC MATRIX INTO AN UPPER TRIANGULAR MATRIX. 


A: 


mM: 
NA: 


500 


10 


600 


20 


503 


22 


3:0 
40 
45 


50 
200 


ARRAY ORIGINALLY CONTAINING MATRIX TO BE DECOMPOSED. 
AT THE END IT CONTAINS THE UPPER TRIANGULAR MATRIX. 


ORDER OF A 
ROW AND COLUMN DIMENSION OF A 
INTEGER N 


DIMENSION A(40,40) 
15 Web) ky doe 
WRITE (IOT, 2) 


FORMAT (’ ZERO OR NEGATIVE RADICAND’ ) 
WRIidE (161,500) Ayia) 
FORMAT (’ Atl, 1) 7 7 Fisers) 
GSO 1OeZ2ce 


A(1,1)=SORT(A(1,1)) 

DO 10 J=2,N 

Reese ey Fal) 

DO 40 I=2,N 
WRITE(IOT,600) I 
FORMAT (13) 

T1i=I-1 

D=A (1,1) 

DO 20 L=1,I11 

D=D-A(L,I1)*A(L,I) 
WREGE (lOr 500) 1,.D Ati 
FORMAT (13,2F13.11) 

ERG ho ele a) le eco 

A(I,1)=SQRT(D) 
IF(I .EQ. N) THEN 

GO TO 45 

END IF 

P2=142 

DO 40 J=12,N 

D=At1 dd) 

DO 30 L=1,I11 

D=D-A(L,I1) *A(L,J) 
WRITE (TOT, 503): 1,D°A(a7 2.) 

AAT di) = D/A 1) 

DO 50 I=2,N 

T1l=I-1 

DO 50 J=1,I11 

A({I,J) =0 

RETURN 

END 


38 


APPENDIX B 


Case Run One 


Branch Flow Rate Pin Pout Re 
(ft/sec) (psf) (psf) 
1 .900035966 .03996 BOS 16:0 1299. 
2 7OOO19242 .05160 ZOLIOS 695. 
3 .00015641 ~-O1905 HOO 74 I 5G oe 
4 .00019242 .00741 .03996 6952 
5 .00016724 ZOOS Ss .03996 604. 
6 .00016724 FOS 160 sO 19 604. 
a 00003601 soLgSOS Home 4 ES) O)e, 
8 .00003601 O07 al .00000 ron 
9 sOOC13 443 ,OO5Ss5 .01719 485. 
EO MOOI 2.5. .01719 .01164 118. 
11 -OO0068S2Z .01164 .00000 249. 
2 .00003281 .00000 PUe S55 Lee 
Case Run Two 
Branch Flow Rate Pin Pout Re 
(m/sec) (N/sm) (N/sm) 

al: BOOO1L3S 296 =-.28579 .37451 698. 
2 MOOOOG6 95 .37451 ,13925 362. 
3 .00005561 PGS eS nOSO52 292. 
4 .00006895 = OS OS Zz aioe? 9 S62 
5 .00006401 - .04323 S265 79 S336u 
6 .00006401 73 (ao .13195 Booe 
q .00001333 213925 BOG y 2 1Oe 
8 HOOO001 333 =f Us) s)% .00000 Oe 
9 .00005134 ~.04323 eS 270. 
0) [OOO001267 st3to5 mm Ohotsie a ove 
11 "OO002600 POSS 72 “O00 0'0 ee 
ily mOOOO126 7 SOV ELONS 1G. .04323 67. 
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